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Abstract 

The generalizations of instantaneous frequency and instantaneous bandwidth to a bivariate signal are derived. These are 
uniquely defined whether the signal is represented as a pair of real-valued signals, or as one analytic and one anti-analytic signal. 
A nonstationary but oscillatory bivariate signal has a natural representation as an ellipse whose properties evolve in time, and 
this representation provides a simple geometric interpretation for the bivariate instantaneous moments. The bivariate bandwidth 
is shown to consists of three terms measuring the degree of instability of the time-varying ellipse: amplitude modulation with 
fixed eccentricity, eccentricity modulation, and orientation modulation or precession. A application to the analysis of data from a 
■ free-drifting oceanographic float is presented and discussed. 

Index Terms 

Amplitude and Frequency Modulated Signal, Analytic Signal, Instantaneous Frequency, Instantaneous Bandwidth, Multivariate 
Time Series. 

^ ; I. Introduction 

THE representation of a nonstationary real-valued signal as an amplitude- and frequency-modulated oscillation has proven 
to be a powerful tool for univariate signal analysis. With the analytic signal [1] as the foundation, the instantaneous 
m ' frequency [2], [3] and instantaneous bandwidth [4] are time-varying quantities which measure, respectively, the local frequency 
, content of the signal and its lowest-order local departure from a uniform sinusoidal oscillation. The instantaneous quantities are 
intimately connected with the first two moments of the spectrum of the analytic signal, and indeed can be shown to decompose 
the corresponding frequency-domain moments across time. 

In the past few years there has been a great deal of interest in nonstationary bivariate signals. A variety of methods have 
been proposed for their analysis, including: examining pairs of time-frequency representations [5], [6], i.e. the dual frequency 
,— 1 ' spectrum, the Rihaczek distribution, and relatives; creating suitable coherences based on these objects [7]; a bivariate version of 
^ , the empirical mode decomposition [8]; local analysis of wavelet time-scale ellipse properties [9]; and an extension of wavelet 
! ' ridge analysis for modulated oscillatory signals [10], [11] to the bivariate case [12]. The statistical and kinematic properties 
, of stationary bivariate signals have also been recently revisited in a number of studies [13]— [16]. 

As an example of the importance of bivariate nonstationary signals, Fig. QJ shows the trajectory of a freely drifting 
, oceanographic instrument [17] as it follows the ocean circulation for hundreds of kilometers. Acoustically-tracked subsurface 
• "floats" such as this one [18] are a prominent data source for studying the physics of the turbulent ocean. There exist hundreds 
ON . of such records from all over the world, and their analysis constitutes an important and active branch of oceanographic research. 

' Panels (b) and (c) in this figure show the estimated oscillatory signal and background using the wavelet ridge method of [12], 
J> , The estimated signal, a modulated bivariate oscillation, is represented in Fig. QJ) as a sequence of ellipses plotted at different 
' times. The analysis of this data will illustrate the great utility of extending the concept of instantaneous moments to the bivariate 
V s * case. 

The goal of this paper is to extend the main tools for analyzing univariate modulated oscillations — the instantaneous 
moments — to the bivariate case. A key ingredient for the analysis of a multivariate signal is the specification of a suitable 
structural model, the purpose of which is to condense the information from a set of disparate signals into a smaller number of 
relevant and intuitive parameters. Work thirty years ago in the oceanographic literature [19]— [23] demonstrated the utility of 
considering a stationary bivariate signal to be composed of a random ellipse at each frequency. Here the recently introduced 
notion of a modulated elliptical signal [12] will be used as the foundation for our analysis. 

The structure of the paper is as follows. The univariate instantaneous moments are reviewed in Section IH-AI Section [III] 
introduces the modulated elliptical signal, and expressions for pairs of instantaneous moments are derived. This suggests, for a 
multivariate signal, defining a joint instantaneous frequency and a joint instantaneous bandwidth which integrate, respectively, 
to the frequency and bandwidth of the average of the individual Fourier spectra. The joint instantaneous moments are derived 
for general multivariate signals in Section |IV] and the bivariate case is then examined in detail. These quantities are shown to 
be independent of unitary transformations on the analytic signals, including in particular real rotations of the coordinate axes. 
An application to the oceanographic data of Fig. Q~|is presented in Section [V] and the paper concludes with a discussion. 
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Fig. 1, As an example of a bivariate time series, panel (a) shows the trajectory of a freely drifting oceanographic float as it follows the ocean currents in 
the northeast subtropical Atlantic. The original data (a) is decomposed into a modulated bivariate oscillation (b) plus a residual (c). In (b), the modulated 
oscillation is represented by "snapshots" at different times, plotted following the instrument location, and black and gray ellipses are alternated for clarity. 
The time interval between snapshots varies in time, and is equal to twice the estimated local period of oscillation. The beginning of the record is marked by 
a triangle. 



The most important results of this paper concern the instantaneous bandwidth, a fundamental quantity which accounts, 
together with variations of the instantaneous frequency, for the spread of the Fourier spectrum about its mean [4]. The correct 
generalization of the univariate bandwidth to the multivariate case has a surprising but intuitive form. In particular, the bivariate 
bandwidth admits an elegant geometric interpretation as a fundamental measure of the degree of instability of a time-varying 
ellipse. It consists of three terms: root-mean-square amplitude modulation, eccentricity modulation or distortion, and orientation 
modulation or precession, all of which contribute to the Fourier bandwidth. 

All data, numerical algorithms, and scripts for analysis and figure generation, are distributed as a part of a freely available 
package of Matlab functions. This package, called JLAB, is available at the first author's website, http://www.jmlilly.net 

II. Background 

This section gives a compact review of the theory of univariate instantaneous moments, drawing in particular from [4], and 
gives definitions which will be used throughout the paper. 

A. The Analytic Signal 

A powerful model for a real-valued nonstationary signal x(t), assumed deterministic and square-integrable herein, is the 
modulated oscillation 

x{t) = a x (t) cos <j) x (t) (1) 

where a x (t) > and <j> x (t) are here defined to be particular unique functions called the instantaneous canonical amplitude 
and canonical phase, respectively [24]. It should be recognized that, in general, the choice of an amplitude and phase pair of 
functions in (T} for a given signal x(t) is not unique [see e.g. 25]. However, the canonical pair is of fundamental importance 
because it separates the signal variability into a x (t) and cos <fi x (t), so that for suitable classes of signals these concur with a 
specified amplitude and phase (see Bedrosian's theorem). 

The canonical pair [a x (t), (f> x (t)} is defined in terms of the analytic signal x + (t), defined as [3], [26] 

1 f°° x(u) 

x+(t) = 2A[x](t) = x(t) +i~ f -^—t-du (2) 

where "£" is the Cauchy principal value integral and A is called the analytic operator. This creates a unique complex-valued 
object x + (t) from any real-valued signal x(t), and permits x(t) to be associated with a unique amplitude and phase through 

a x (ty"" lt) =x+(t), x+(t)^Q (3) 

which defines the canonical pair [a x (t), 4> x (t)]. The original signal is recovered from the analytic signal via x(t) — 5Ji {x+(t)}, 
where "5R" denotes the real part. 
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The action of the analytic operator is more clear in the frequency domain. X(ui) is the Fourier transform of x{i) 

x(t) = — / X{u)e wt duj (4) 

27T J-oc 

and likewise X+(cu) is the Fourier transform of x+(t). The time-domain operator in <(2j becomes in the frequency domain 
simply [27] 

X + (lu) = 2U{uj)X{uo) (5) 

where U(u>) is the unit step function. Thus in the frequency domain, the analytic signal is formed by doubling the amplitudes 
of the Fourier coefficients at all positive frequencies, while setting those of all negative frequencies to zero. 

B. Global and Instantaneous Moments 

Using the analytic signal, one may describe both the local and the global behavior of any real-valued signal x(t) as a 
modulated oscillation. A global description is afforded by the frequency-domain moments of the spectrum of the analytic 
signal. The global mean frequency and global second central moment of x(t) are defined, respectively, by 

1 



^ u\X + {uj)\ 2 duj (6) 

ATtCx Jo 
1 f°° 

w * = I (w-w:0 2 pr + (w)| 2 dw (?) 



where 



-i pOO />oo 

— J IX+HI 2 dcj = J \x + (t)\ 2 dt 



= — / \X+(u)\'du= I \x+(t)\*dt (8) 



is total energy of the analytic part of the signal. The second central moment ~a x characterizes the spread of the spectrum about 
the global mean frequency; thus a x is also known as the bandwidth. 

It would be greatly informative to relate the global moments oJ x and a x to the time evolution of the signal. This leads to 
the notion of instantaneous moments, time-varying quantities which recover the global moments through integrals of the form 
[4], [28] 

/OO 
al(t)u> x (t)dt (9) 
-oo 

*l = S- 1 / a 2 x (t)a 2 x (t)dt. (10) 



It is a fundamental result that the derivative of the canonical phase 

u x {t) = j t <t> x {t)=%^j t \nx + {t)X (11) 

satisfies (0, and is termed therefore the instantaneous frequency [2], [3]. It is important to point out that the instantaneous 
moments u x (t) and er^(£) are not themselves moments; rather, they are instantaneous contributions to global moments. 

The instantaneous second central moment cr^ (t) is not uniquely defined by the integral (fTOl . because more than one function 
of time may be shown to integrate to the corresponding global moment. However, the constraint that a x (t) be nonnegative 
definite, like W 2 ,, leads to the unique definition [4] 



| j- t x+(t) - iu x x + (t) 



°U*J = i 77775 ■ ( 12 ) 



The square root of the instantaneous second central moment, a x (t), then has an intuitive interpretation as the average spread 
of the frequency content of the signal about the global mean frequency uJ x . 
It has been shown by [4] that the fractional rate of amplitude modulation 

^ (t)s dlB^W =R fd lnMt) l (13) 



dt [dt 

plays an important role within the instantaneous second central moment. After a simplification, (fT2l becomes 

a 2 x (t) = [u; x (t)-uJ x } 2 + v 2 x (t) (14) 

from which it is clear that variations of the instantaneous frequency uo x (t) about the global mean frequency ZJ X do not account 
for the entirety of the global second central moment a x . The remainder is accounted for by v x (t), and since a x is known as 
the global bandwidth, v x (t) is called the instantaneous bandwidth [4]. 
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Sketch of ellipse 
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Fig. 2. A sketch of the parameters in the modulated elliptical signal, as described in the text. 

III. The Modulated Elliptical Signal 
This section introduces a model for a modulated bivariate oscillation, extending and simplifying the development in [12]. 

A. Modulated Ellipse Model 

Bivariate signals will be represented both in the complex-valued form z(t) = x(t) + iy(t), and in vector form x(t) = 
[x(t) y(t)~\ , with x(t) and y(t) both real-valued and where "T" denotes the matrix transpose. Our starting point for 
describing bivariate oscillations will be the modulated elliptical signal [12] 

z(t) = e l8{t) {a(t) cos <f>(t) + ib(t) sin <j)(t)} (15) 

which represents z(t) as the position traced out on the [x,y] plane by a hypothetical particle orbiting a time-varying ellipse; 
see the sketch in Fig. [2] The two angles 9(t) and <fr(t) are defined on the principal interval (— n, ir), while a(t) > 0; b(t) may 
take on either sign for reasons to be see shortly. 

The ellipse has instantaneous semi-major and semi-minor axes a(t) and \b(t)\, and an instantaneous orientation of the major 
axis with respect to the x-axis given by 8(t). The angle <fi(t), called the orbital phase, determines the instantaneous position 
of the particle around the ellipse perimeter with respect to the major axis. While a(t) is nonnegative, the sign of b(t) is 
chosen to reflect the direction of circulation around the ellipse. Note that the modulated univariate signal (fl} is included as the 
special case b(t) — 0, 9(t) — 0. Thus we can view the modulated elliptical signal as an extension of the amplitude/frequency 
modulated signal to the bivariate case. 

The meaning of the "instantaneous" ellipse properties referred to in the previous paragraph is more clear if we separate the 
ellipse phase from the other variables. Introduce 

z(t, t') = e l9[t) {a{t) cos [0(t) + u)^{t)t'] + ib(t) sin [<j){t) + (*)*']} (16) 

where t' is seen as a "local" time. With t fixed, z(t, t') continually traces out the perimeter of a "frozen" ellipse as t' varies; the 
period of orbiting the ellipse is 2Tr/oj^(t). It is these frozen ellipses which are plotted in Fig. [TJi for the estimated modulated 
elliptical signal corresponding to the time series in Fig. 0J. Note that z(t) itself will not in general trace out an ellipse if the 
ellipse geometry varies in time. However, if the parameters of the ellipse geometry a(t), b(t), and 6(t) are slowly varying with 
respect to the phase 4>(t), then z(t) will approximate an ellipse. 

B. Analytic Signal Pairs 

The parameters of the modulated elliptical signal are made unique for a given pair of real-valued signals by referring to a 
pair of analytic signals, just as the univariate modulated signal model is made unique by referring to a single analytic signal. 
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In vector notation the modulated elliptical signal ( fT3T > becomes 

\x(ty 



j[0(t)]n 



i<t>{t) 



a(t) 
-ib(t) 



where 



3(6) 



cos f 
sin t 



— sm( 
cos 8 



is the counterclockwise rotation matrix. Thus we may write 

e^ (t) J [6{t)] 







x+(t) 




~2A[x](t) 


-ib(t) 








2A[y](t)_ 



(17) 



(18) 



(19) 



and this equation is taken as our definition of the ellipse parameters in terms of the analytic portions of x(t) and y(t). In terms 
of the analytic signals x + (t) and y+(t), 



z(t) = $l{x + (t)} + m{y + (t)} 



(20) 



recovers the original complex-valued signal. 

We emphasize that a(t), b(t), 9(t), and 4>(t) appearing in ( fl5] l and ( fT71 ) are not arbitrary functions which one is free to 
specify. They are, instead, unique properties of a given bivariate signal. While any bivariate signal can be represented as a 
modulated ellipse, just as any univariate signal can be represented as a modulated oscillation via the analytic signal, it is not 
always sensible to do so. If one begins with a pair of real-valued signals x(t) and y(t) which are, say, finite samples of a 
noise process, the ellipse parameters would likely not yield any sensible information; this would inform you that the modulated 
ellipse representation is not appropriate. 

One may also associate a second pair of analytic signals with z(t), namely 

(21) 

the components of which are the analytic parts of z(t) and of its complex conjugate. Then 

z(t) = z+(t) + z*_(t) (22) 

decomposes z(t) into counterclockwise and clockwise rotating contributions, respectivelyQ Note that for convenience we have 
defined to be an analytic, rather than an anti-analytic, signal. Making use of the unitary matrix 

T = -^= 

one finds 



" z+(t) ' 




~A[z](t)- 


z_{t) 




A[z*](t)_ 



(23) 



z+(t) 

Z.(t) 



1 



V2' 



x+(t) 

y+(t) 



(24) 



as the relationship between the two different pairs of analytic signals. In order to distinguish between these two pairs, we refer 
to x+(t) and y+(t) as the Cartesian pair of analytic signals, and to z + (t) and z_(t) as the rotary pair (borrowing a term 
from the oceanographic community [19], [20]). 

The ellipse parameters have simple expressions in terms of the amplitude and phases of the rotary analytic signals. Letting 



z+(t) = o+^e^+W 
z-(t) = a_(t)e^- (t) 



(25) 
(26) 



and comparing the rotary decomposition ( 1221 with T times the matrix form of the modulated ellipse model ( fT9] l, one finds 
[12] 



a(t) = a+(t)+a-(i) 

b(t) = a+(t)-a-(t) 

Cj>{t) EE [0 + (t) + ^_(t)]/2 

0(t) EE fo + (t)-^_(t)]/2 



(27) 
(28) 
(29) 
(30) 



as expressions for the ellipse parameters in terms of the amplitude and phases of the rotary analytic signals. The ellipse 
parameters are thereby uniquely defined. Unique relations of the ellipse parameters to the parameters of the Cartesian analytic 
signals are then implied by d24"i >. and are given in Appendix [A] 

'The factor of two difference in defining analytic versions of real- valued and complex- valued signals, i.e. the last equality in )19t versus j2U . prevents 
factors of two from appearing in the decomposition equations )20t and i22i . 
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C. Amplitude and Eccentricity 

It is convenient to replace the semi-major and semi-minor axes a(t) and b(t) with the root-mean-square amplitude 



'a 2 (t) + 


b 2 (t) 








b 2 (t) 


a 2 (t) + 


b 2 (t) 



«(*) = \l w „ - (31) 
together with a relative of the eccentricity [29], [12] 

A( * )sr *5# + ^ (32) 

which measures the signal's departure from circularit)0. Here the direction of rotation of the "particle" around the time-varying 
ellipse is denoted by 

r z = sgn{b(t)} (33) 

which we assume henceforth to be constant over the time interval of interest. These definitions give 



a(t) = K (t)y/l + \\{t)\ (34) 
b(t) = r zK {t)y/l-\X(t)\ (35) 

as expressions for the semi-major and (signed) semi-minor axes. Note that r z is well defined for purely circular signals for 
which a(t) = \b(t)\, whereas sgn{A(t)} vanishes. 

The magnitude of A(t), like the eccentricity ecc(t) = v/l — b 2 (t) / a 2 (t) , varies between |A(t)| = for purely circularly 
polarized motion and |A(i)| = 1 for purely linearly polarized motion. We refer to |A(i)| as the ellipse linearity and to X(t) as 
the signed linearity. While not in common use as a measure of eccentricity, we prefer X(t) to other such measures since it 
leads to simple forms for subsequent expressions. 



D. Rates of Change 

The rates of change of the ellipse parameters have special interpretations. The time derivative of the orbital phase, termed 
the orbital frequency Lj$(t) = 4l 4>{t), gives the rate at which the particle orbits the ellipse, while the time derivative of the 
orientation angle is the precession rate u>e(t) = 4? 8(t). Amplitude modulation of the ellipse involves a time derivative of k(£), 
while variation in X(t) means that the shape of the ellipse is distorting with time. However, the ellipse parameters and their 
derivatives do not have an immediately evident relationship to the global moments of z(t). Thus, while these are useful local 
descriptions of joint signal variability, they are not interpretable as instantaneous moments. 



E. Instantaneous Moment Pairs 

In Section ITl-AI it was shown that the instantaneous moments of a univariate signal — instantaneous frequency and bandwidth, 
in particular — provide a powerful description of local signal variability with a direct relationship to the global signal moments. 
We aim to identify the analogous quantities for a bivariate signal. 

A bivariate oscillatory signal can be equivalently expressed in terms of two pairs of analytic signals, the Cartesian pair 
x + (t) and y+(t) or the rotary pair z+(t) and z~(t). Each of these, in turn, can be analyzed by their individual instantaneous 
moments, using the ideas described in Section III-AI In particular, these two pairs of analytic signals lead to two pairs of 
instantaneous frequencies [u! x (t), u> y (t)] and [u>+(i), and two pairs of instantaneous bandwidths [v x (t), v y (t)] and 

[v+(i), v-(t)], which contribute to two corresponding pairs of global moments. This approach is, however, unsatisfactory. 
Rather than a unified description of signal variability, an analysis based on these pairs of instantaneous moments achieves a 
description of two disparate halves of the bivariate signal z[t) considered separately. 

Between the Cartesian and rotary analytic signals, one can represent the signal as a pair of oscillations in two different ways. 
The difficulty of a pairwise decomposition is highlighted by considering a range of different values of the linearity |A(£)|. 
Generally speaking, the decomposition of z(t) into the rotary pair of analytic signals is more appropriate when one of z + (t) 
and Z-{t) is much larger than then other, e.g. |.s+(t)| S> On the other hand, when (in some rotated coordinate system) 

one has |x+(i)| 3> |y+(i)|, then the Cartesian decomposition is more appropriate. These two cases correspond to |A(f)| ~ 
and |A(t)| ~ 1, respectively. In between these two cases, or for a signal ranging over both extremes, there is a gray area in 
which neither decomposition is particularly appropriate. A unified description is therefore necessary. 

2 Note that herein we always use the term "circular" to describe the shape traced out on the [x,y] plane by a bivariate signal, rather than in the statistical 
sense of a circularly symmetric or proper complex-valued signal [e.g. 13], 
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To begin with, the instantaneous moments of the analytic signal pairs can be cast in terms of the ellipse parameters, to show 
how changes in the ellipse geometry are expressed in time variability of the instantaneous moments pairs. Inverting (|27H30l l, 
and utilizing d34l i and d35l l, leads to 

0±(i) = 0(t)±0(f) (36) 

a ± (t) = ^^l± rz ^>?{£) (37) 

giving the parameters of the rotary analytic signals in terms of the ellipse parameters. (In these and subsequent expressions, 
the signs on the right-hand-side are understood to be chosen to match the signs on the left.) Note that the amplitudes satisfy 
+ = K 2 (t), and we may identify K 2 (t) as the sum of the instantaneous power of the two rotary signals. When 

r z > 0, and the ellipse rotates in the positive (counterclockwise) direction, we have a+(t) > o_(t) as expected, whereas the 
opposite is true for negative rotation. Note that as |A(f)| approaches zero, and the signal becomes nearly circular, the smaller 
of o+(t) and ti-(t) approaches zero while the larger approaches n(t). 
Differentiating ( f36l > and the logarithm of (f3Tb leads to 

w±(i) = w*(t)±w e (*) (38) 

as expressions for the rotary instantaneous frequencies and bandwidths. It is useful to note an asymmetry of the rotary 
bandwidths. For small |A(t)| with r z > 0, the denominator of the second term in u_(t) is close to zero while that of v+(t) is 
close to unity; the situation is reversed for r z < 0. This shows that for small |A(t)|, the bandwidth of the weaker of the two 
rotary signals is much more sensitive to variations in the degree of eccentricity than is the bandwidth of the stronger signal. 

As discussed above, these expressions are useful in the near-circular case |A(f)| ~ 1. The same procedure carried out in 
terms of the Cartesian analytic signals x + (t) and y+(t), presented in Appendix lAl leads to more complicated relationships 
between variations of the ellipse geometry and the instantaneous moments of the two real-valued signals x(t) and y(t). In 
particular, we may note in those relationships the explicit dependence on the ellipse orientation 6(t). This renders the Cartesian 
instantaneous moments quite problematic for a unified description of the signal variability, since a simple coordinate rotation 
will cause these quantities to change. 

There is therefore a disconnect between the moment-based description of bivariate signal variability, grounded on the analytic 
signals, and the modulated ellipse model of joint structure. This motivates the development of the next section. 

IV. Joint Instantaneous Moments 

The definitions of instantaneous moments in Section Hi- Al which are standard, can be extended to accommodate multivariate 
oscillatory signals. Let 

x+(i) = [a? +i i(t) x + , 2 (t) ... x + . N {t)f (40) 

be a vector of N analytic signals, and let X+(cj) be the corresponding frequency-domain vector; these are related by the 
inverse Fourier transform _ 

x+(t) = — y e^X+H^. (41) 

Here we begin with the general multivariate case and subsequently examine the bivariate case, N = 2, in detail. 
A. Joint Moments 

It is presumed that the components of x+(i) are sufficiently closely related that we should seek a unified description of 
their variability. To this end it is reasonable to define the joint analytic spectrum as the normalized average of the spectra of 
its N components, 

= OlX+HH 2 (42) 



where ||x|| = Vx^x is the Euclidean norm of a vector x, "i?" indicating the conjugate transpose, and where 

5 X =— / ||X + (w)|| 2 dw= / \\x+(t)fdt (43) 
2tt Jo J-oo 

is the total energy of the multivariate analytic signal. One may then define the joint global mean frequency 

' ' cuS x (uj)duj (44) 



27r£ x jo 
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oZ = -±- I {uj-uj^f S x (co)dco (45) 



associated with the joint analytic spectrum, together with 

1 

2tt£ x jo 

which is the joint global second central moment. 

The joint instantaneous frequency cu x (t) and joint instantaneous second central moment <J x (t) of x + (t) are then some 
quantities that decompose the corresponding global moments across time, i.e. which satisfy 

/oo 
\\ X+ (t)\\ 2 uj x (t) dt (46) 
-OO 

/OO 
\\ X+ (t)\\ 2 <jl(t)dt (47) 
-oo 

noting that ||x + (i)|| 2 is the instantaneous total analytic signal power. Since the univariate instantaneous frequency ( fTTT i may 
be rewritten 

^ (t)59 (jj£M)\ = 9 H (t)f y 1 } m 

l * J k+(t)l 

we define the joint instantaneous frequency according to the same form, 

9 {*?(*)*#} 

and find that this does indeed satisfy (l46b . 

Likewise, the second instantaneous central moment, which we define as 

2 m _ ||| x +W-^xX + (Q|| 2 

ff x' = ii 7TVM2 (5U) 

is a nonnegative-definite quantity satisfying (1471 . This gives the normalized departure of the rate of change of the vector-valued 
signal from a uniform complex rotation at the constant single frequency ZU X . Equations (|49l i and ( T50b are clearly the natural 
generalizations of (fTTT > and ([T2l to multivariate analytic signals. 



B. Joint Instantaneous Bandwidth 

Furthermore we may generalize the notion of instantaneous bandwidth to a multivariate signal. We define the joint instan 
taneous bandwidth via its relationship to the joint second central moment 

i2 



<{t) 



<7 x (t) - - LO y 



(51) 



by comparison with (TBI for the univariate case. That is, we define the squared instantaneous bandwidth to be that part of the 
instantaneous second central moment not accounted for by deviations of the instantaneous frequency from the global mean 
frequency. 

This definition of the bandwidth leads, after some manipulation, to 

,2 



[^x+(t)-«^(t)x+ (t)| 
l|x+(t)|| 2 



(52) 



which is the normalized departure of the rate of change of the vector-valued signal from a uniform complex rotation at a single 
time-varying frequency ui x (t). For x+(f) a 1 -vector consisting of a single analytic signal, x+(t), (l52l becomes 

i2 



[(*) = 



— In a^i) 
at 



= v%{t) 



(53) 



and the joint instantaneous bandwidth correctly reduces to the univariate bandwidth defined in (fLTt . 

The contributions to the instantaneous second central moment and the instantaneous bandwidth are perhaps more clear if 
we write (f5Qb and (l52l out as summations 



^(*) 



Sn=l°n(*){^W + KW-57x| 2 } 



Etl ffl nW 



Eti«iM{^W + k(f)-^)| 2 } 



En=l a n(*) 



(54) 

(55) 
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where a n (t) is the amplitude of the nth analytic signal, and so forth. The first of these states that amplitude modulation, as 
well as departures of the instantaneous frequencies from the global mean frequency ZU X , contribute to the second instantaneous 
central moment. The second states that amplitude modulations together with departures of the instantaneous frequencies from 
the time-varying joint instantaneous frequency oj x (t) contribute to the squared joint instantaneous bandwidth. In both cases 
contributions from different times and different signal components are weighted according to the instantaneous power a^(t). 



C. Invariance 

The joint instantaneous moments defined above have the important property that they are invariant to transformations of the 
form 

y+(t) = cUx+(i) (56) 

where c is some constant and U is an N x N unitary matrix; these may be termed scaled unitary transformations. Since in 
the frequency domain we have also Y + (uj) = cUX + (cj), it is obvious that 

YffMY + M 

= ^J-m Y +( w ) Y +H^ =S *(w) (57) 

hence the joint analytic spectrum is unchanged, as are the joint global moments ZJ y = ZJ X and ~o 2 y = er 2 . The joint instantaneous 
frequency transforms as 

" y (t) = b7W " x ? (t)x +W - " X (t) (58) 

and hence remains unchanged. Similarly it is easy to see that the joint instantaneous bandwidth u x (i) and joint instantaneous 
second central moments cr^(t) and are also all invariant to scaled unitary transformations. In particular, these joint instantaneous 
moments are invariant to coordinate rotations. 



D. Joint Bivariate Moments 

We now examine the joint instantaneous moments for bivariate signals by taking either x+(t) = [.z+(i) z_(t)l or x+(t) = 
y+(t)\. Since these two vectors are related by the scaled unitary transformation (|24t . from the invariance shown in 
the preceding subsection, the results will be identical in either case. The bivariate joint instantaneous moments integrate to the 
global moments of the joint analytic spectrum 

^ |x + H| 2 + |y + H| 2 = |z + H| 2 + |z_M| 2 (59) 

£ x + £y £ . + £ 

where £ x = ^ J 00 |X + (o;)| 2 da;, and so forth. Note that we will use the subscript "z" rather than the subscript "x" to 
specifically denote bivariate quantities. 

To see d59l >. it is helpful to relate the joint analytic spectrum to the two-sided spectrum of the complex-valued signal z(t), 
which has a Fourier transform Z(lo). The Fourier transforms of z + (t) and z_(t) are, respectively, 

Z+(lu) = U{u)Z{u) = ^[X + (u)+iY+{u)] (60) 

Z_(w) = U{u>)Z{-u) = ^[X + {u)-iY + {u)\ (61) 

where Z(uS) is the Fourier transform of z(t) and where U(uo) is again the unit step function. Inserting these into the second 
and third expressions in (|59l l, we note at once a cancelation of cross-terms, and the equality follows. Note that d60T - l6lT l show 
that z + (t) and have Fourier coefficients drawn entirely from the positive-frequency and negative-frequency sides of 

Z(lu), respectively (hence the notation "+" and "-")■ The joint analytic spectrum simply consists of averaging the positive and 
negative frequency halves of |Z(w)| 2 , followed by a normalization to unit energy. 

1 ) Instantaneous Frequency: To obtain the bivariate instantaneous frequency, we begin with the definition d49l i and insert 
expressions ((37) and ( T38T > for the rotary instantaneous frequencies and amplitudes. This leads to 



w z (t) = (t) + r z v/l-A 2 (t) cog (t) (62) 

for the bivariate instantaneous frequency written in terms of the ellipse parameters. Alternatively, we could have used expressions 
(IToTfTTT i and d74l - [75l ) from Appendix |A] for the Cartesian analytic signals. Inserting these into (|49| | we obtain again d62l l. as we 
must on account of the invariance proved in the last subsection. The relative difficulty of directly manipulating the cumbersome 
Cartesian expressions illustrates the importance of this general result. 

The joint bivariate instantaneous frequency uj z (t) gives a correct measure the time-varying frequency content of a bivariate 
signal regardless of the polarization state. For |A(t)| = 0, the signal is purely circular, and ui z (t) becomes u>^{t) + r z uig(t). 
Comparison with (l38l l shows that this is u>+(t) if r z > and uJ-(t) if r z < 0, that is, for circularly polarized motion uo z (t) 
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Fig. 3. An ellipse with uniformly increasing relative amplitude (a), i.e. v K (t) /uj z (t) constant, an ellipse with uniformly increasing eccentricity (b), i.e. 
v\ (t ) /u z (t ) constant, and an a uniformly precessing ellipse (c) with vg(t)/u> z (t) constant. The three constants have been chosen to all be the same value 
of 0.025. The bold portion of the line in all three panels show an initial single orbit. A circle marks the beginning of each record and an "x" marks the end 



becomes the instantaneous frequency of the non-vanishing rotary component. On the other hand, in Appendix [A] it is shown 
that for a purely linear signal, |A(i)| = 1 or vanishing minor axis \b(t)\ = 0, the Cartesian instantaneous frequency along the 
coordinate axis aligned with the ellipse major axis is equal to u)^,(t), but (f62t shows that with |A(t)| = 1 then this is also equal 
to Ld z (t). Thus the use of u> z (t) is quite desirable as it appropriate for any polarization, unlike the partitioning into a rotary 
pair or a Cartesian pair. 

2) Instantaneous Bandwidth: Similarly the bivariate bandwidth, expressed in terms of the ellipse parameters, is 



v 2 At) 



d\i\K,(t) 



dt 



1 dX(t) 

2 dt 



A 2 (*H 2 (i) 



(63) 



1-A 2 (i) 

as follows from the definition d53l together with (f39l > or (|79H80l l. This expression is quite illuminating. The first term in 
measures the strength of amplitude modulation, the second is the rate of ellipse distortion on account of changing eccentricity, 
and the third is due to the precession of the ellipse. Note each of the these quantities is independent of coordinate rotation, 
that is, the orientation angle 6(t) does not explicitly appear. We define these three terms as 



v 2 z (t) = v 2 K (t) + v 2 x (t) + vl(t) 



(64) 



which we call the squared amplitude bandwidth, deformation bandwidth, and precession bandwidth, respectively. 

An illustration of the bivariate bandwidth is presented in Figure [3] Three different time-varying elliptical signals are plotted, 
all having constant orbital frequency ui$(t) and each with exactly one of the three terms in (l63l being nonzero. The quantities 
v K {t), v\(t), vg(t) are each set equal to the value of 0.025xcu^(t) for the three cases respectively. It is clear that if any of 
these terms were to become too large, the usefulness of our description of the signal — as an ellipse the properties of which 
evolve with time — would become questionable. Thus v z (t) quantifies the lowest-order departure of the bivariate signal from 
periodic motion tracing out the periphery of a fixed ellipse. 

The bivariate bandwidth is therefore a fundamental quantity reflecting the degree of instability of the elliptical motion. It is 
remarkable and surprising that the definition d49| i, which is a power-weighted average of corresponding univariate quantities, 
should be identical with d63l which is clearly an expression of the instability of elliptical motion. That these expressions are 
equivalent underscores the fact that bandwidth is itself a measure of oscillation stability, and conversely, that the degree of 
instability of an elliptical signal is interpretable as a bandwidth. 

The deformation bandwidth v \{t) deserves further comment. Note that this quantity can be expressed in the equivalent forms 



M*)l 



d\{t) 



y/1- A 2 (i) dt 

as may readily be verified. These lead to the following approximations 



1 d 



A(t) dt 



Vl-A 2 (t) 



\vx(t)\ 



d\(t) 



dt 



dt^ 1 



[l + 0(A 2 (i))] 

^\W) [1 + 0(1-A 2 (i))] 



(65) 

(66) 
(67) 
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appropriate for the near-circular case |A(i)| ~ and the near-linear case |A(i)| « 1, respectively. In the near-circular case, the 
deformation bandwidth is due to the (small) departures of the linearity |A(i)| from zero, while in the near-linear case it is due 
to the (small) departures of the linearity from unity. Another interpretation may be found by noting, using 



K(*)l = 



a(t)b(t) 



■In 



a(t) 



(68) 



a 2 (t) + b 2 {t) dt 

which states that the deformation bandwidth is due to fractional changes in the ellipse aspect ratio b(t)/a(t), weighted in 
proportion to the ratio of the ellipse area Tra(t)\b(t)\ to the root-mean-square radius. 



V. Application 

This section presents an application of the bivariate instantaneous frequency and bandwidth to a typical signal from physical 
oceanography. 



A. Data and Method 

As discussed already in the introduction, Fig. QJ presents data from a subsurface, acoustically-tracked oceanographic float 
[18]. The data was downloaded from the World Ocean Circulation Experiment Subsurface Float Data Assembly Center 
(WFDAC) at http://wfdac.whoi.edu This particular record, from an experiment in the eastern North Atlantic Ocean that 
is well-known among oceanographers [17], [30], [31], was recorded somewhat to the south and west of the Canary Islands. 
The float was drifting at a depth between 1000 m and 1300 m. Its horizontal position was recorded once per day, as inferred 
by triangulation of acoustic travel times from several nearby sound sources. The region shown ranges from 20-25°N and 
27-22.5°W. 

The loops in the trajectory are known to be the imprint of one of the large "eddies" or "vortices" [32] that are common in 
the ocean, a distant relative of the swirls one observes when stirring coffee; note that the loops seen in Fig. QJ are up to fifty 
kilometers in diameter. An instrument trapped in such a vortex will record both nearly circular motion from orbiting around 
the vortex center, which is well modeled as a modulated elliptical signal as shown by [12], as well as translational or advected 
motion of the vortex center itself. This suggests modeling the observed signal z^{t) = x^(t) + iy^(t) via the unobserved 
components model [33] 

z {o} (t) = z(t) + z {e} {t) (69) 

where z{t) is a modulated elliptical signal and where z^ e '(t) is a residual defined to include everything else. The residual 
z^(t) is expected to include the turbulent background flow together some measurement noise, but may be considered "noise" 
at present since we are interested only in z(t). 

A means of estimating z(t) using an extension of wavelet ridge analysis [10], [11] was developed by [12]. The details of this 
method are not particularly important here, except to note that it is expected to give a reasonable estimate z(i) of the presumed 
(unobserved) modulated elliptical signal z(t). We apply this method, with parameter choices noted in Appendix iBl to the record 
in Fig.QJ. This leads to the estimated bivariate oscillation ~z(t) shown in Fig.QJ), and — by subtraction, z^(t) = z^ (t) — z(t) — 
to the estimated residual signal z^{t) shown in Fig.QJ. In Fig.OJ the signal z(t) is represented by snapshots of the modulated 
signal at different times, by letting the orbital phase vary with the ellipse geometry held fixed, as described earlier in the 
discussion of d!61 l. Note that in Fig. [Tfc) virtually all the looping motions have been removed, leaving behind a large-scale 
meander plus small-scale irregularities. 



B. Instantaneous Amplitude and Frequency 

Having achieved this decomposition, we now analyze in more detail the properties of the estimated bivariate oscillatory 
signal z(t) using the joint instantaneous moments. The estimated signal, shown in Fig. @k together with the instantaneous 
root-mean-square amplitude n(t), exhibits both substantial amplitude as well as frequency modulations; note that we will 
denote all properties of z(t) with a hat, to distinguish them from the properties of the unobserved true signal z(t). In 

particular, during the last third of the record, say after yearday 250, the signal amplitude is greatly reduced compared to the 
earlier time period. This accounts for the transition from large to small ellipses seen in Fig. QJ. 

This signal rotates in a clockwise fashion: the negative rotary component dominates the positive component z+(t), thus 
the signed linearity X(t) (not shown) is everywhere negative; this sense of rotation can be inferred from the phase shift between 
the real and imaginary parts of in Fig. |4^. The linearity \X(t)\ itself (Fig. O?) is generally very small, corresponding to 
nearly circular motion, apart from a few excursions to higher values. Thus only a handful of the ellipses shown in Fig. QJ) 
exhibit substantial eccentricity. 

The instantaneous frequency content of the signal is shown Fig.|4j;. Here, the bivariate instantaneous frequency of u) z (t), 
is presented together with its orbital frequency u)^(t) and precession rate tDg(t). The precession rate is seen to be much smaller 
than the orbital frequency, and consequently, it follows from (l62l that the bivariate instantaneous frequency Q z (t) should 
be close to the orbital frequency uj^t). A transition is seen around around yearday 250 to higher- frequency oscillations, 
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Fig. 4. Instantaneous moment analysis of the estimated modulated elliptical signal z(t), derived from the bivariate time series of Fig. [TJi as described in the 
text. Panel (a) shows the real and imaginary parts of the signal z(t), x(t) (thin solid line) and y(t) (dashed line), together with the associated RMS amplitude 
TZ(t) (heavy solid line). In panel (b), the linearity |A(t)| is shown. Three time-varying frequencies are shown in panel (c), the joint instantaneous frequency 
2 z (t) (heavy solid line), the orbital frequency (thin solid line), and the precession rate 23g(t) (thin dashed line). Panel (d) presents the three terms in 

the bivariate instantaneous bandwidth, the amplitude bandwidth v K (t) (heavy solid line), the deformation bandwidth V\(t) (thin solid line), and the precession 
bandwidth vg(t) (dashed line); in (e), the bandwidths are presented again, but this time divided through by <D z (t)| in order to render them nondimensional. 



corresponding to the transition to smaller amplitudes seen in Fig. |4^. Interestingly, in the higher-frequency portion of the 
record, the precession rate exhibits a tendency for negative rotation, but establishing whether or not this result is statistically 
significant is beyond the scope of the present paper. The precession rate u)g(t) and orbital frequency tD^t) both present some 
jagged or irregular variability, which is likely due to the effect of measurement noise; this interpretation is supported by the 
fact that the irregular variability is particularly present after yearday 250, when the signal strength has weakened. 

Measurement noise is expected to strongly effect tu^(t) and uj$(t) in a case such as this one, when then unobserved signal of 
interest zit) is nearly circular in nature. Here we have ^> and so the observed (noisy) signal z^(t) will have a 

much lower signal-to-noise ratio at positive frequencies than at negative frequencies. The estimated positive rotary signal z+(t) 
is thus expected to be considerably noisier than the negative rotary signal z_ (t), and likewise for the associated instantaneous 
frequencies u)+(t) compared with cD_ (t). Since u)${t) and u>e(t) are, respectively, the sum and difference of uj+(t) and 
noise in uj+(t) will effect them both. 

By contrast Q z (t) is an average over the instantaneous frequencies of two independent signal components. It is therefore 
expected to be more robust against the effects of noise than and cD_(t), and also than the Cartesian analytic signals 

Q x (t) and u) y (t) that would be observed along the coordinate axes. This suggests that although Q<j,(t) and u>g(t) are of great 
interest from the point of view of the modulated ellipse model, lj z (t) provides a superior estimate of the time-varying frequency 
content of a noisy bivariate oscillatory signal. The preceding discussion has also brought up the importance of understanding 
the effects of noise on the estimation procedure, a task which is currently underway. 

C. Instantaneous Bandwidth 

The three terms in the bivariate instantaneous bandwidth of the estimated signal z(t) — v K (t), v\(t), and vg(t) — are shown in 
Fig. |4}l. During most of the record, the amplitude bandwidth v K (t) is larger than the other two. Since these contribute as their 
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squares to the bivariate bandwidth v z (t), we see that in fact the amplitude bandwidth v K (t) accounts for the large majority of 
the estimated joint instantaneous bandwidth v z (t) during most of the record. Thus the lowest-order departure of the bivariate 
signal from a uniform oscillation at frequency u) z (t) is generally due to modulation of the root-mean-square ellipse amplitude, 
with its eccentricity and orientation held fixed. This result could have different important physical interpretations; however, 
such questions are outside the scope of the present paper. 

The three bandwidth quantities are again shown in Fig. [5, but this time each has been divided by uj z (t), thus comparing the 
contributions to the instantaneous bandwidth with the instantaneous frequency. Whereas Fig. [4}i suggests that the instantaneous 
bandwidth increases after yearday 250, Fig.[4£ shows that this transition is offset by the corresponding increase in instantaneous 
frequency. This means that the variability of the modulated elliptical, on time scales proportional to the local period of oscillation, 
is in fact similar during the different periods of the record. If anything, the "relative bandwidth" is somewhat lower during the 
last third of the record then in the remainder. 

VI. Discussion 

This paper has examined the extension of instantaneous moments to bivariate signals. To accomplish this, "joint instantaneous 
moments" were defined for a general multivariate signal. These are time-varying quantities that integrate to the global moments 
of the average analytic spectrum, that is, the spectrum of the analytic part of the signal components averaged over the number 
of components. While the joint instantaneous frequency is simply the power-weighted average of the instantaneous frequencies 
of the components, the joint instantaneous bandwidth has an unexpected but very intuitive form. It measures the extent to 
which a multivariate oscillation does not evolve simply by oscillating at the single, time varying frequency given by the joint 
instantaneous frequency. 

These joint instantaneous moments together with the notion of a modulated elliptical signal represent a powerful means for 
analyzing bivariate signals. The "modulated elliptical signal" is a representation of a bivariate oscillation, alternately regarded 
as a pair of analytic signals, as a single time-varying structure. The joint moments describe the nature of the variability of this 
structure, combining the partial descriptions achieved from the individual moments of the members of a signal pair. Expressing 
the bivariate instantaneous frequency and bandwidth in terms of the time-varying ellipse parameters leads to illuminating 
expressions. In particular, the instantaneous bandwidth was shown to consist of three terms: amplitude modulation, deformation 
or eccentricity modulation, and orientation modulation or precession. These three quantities thus express the three basic ways 
ellipse geometry can change with time; the bivariate bandwidth is therefore seen as the degree of instability of a modulated 
elliptical signal. 

An application of the instantaneous moments to a bivariate time series from oceanography leads to a wealth of information 
regarding the time variability of the signal, which in this case could reflect time variations in an underlying oceanic vortex 
structure. In addition to pursing physical questions raised here, the most important outstanding tasks involve quantifying the 
errors involved in estimating bivariate oscillatory signals. There are both "deterministic" errors due to the time variability of 
the signal, which have recently been examined for univariate wavelet ridge analysis [34], as well as random errors due to 
background noise. These topics are the subjects of ongoing research. 

Appendix A 
The Cartesian Analytic Signals 
Here we derive expressions for the Cartesian instantaneous moments, that is, the instantaneous moments of x + (t) = 
o-x (t) e 1 * 1 and y+(t) = a y (t) e^y^ for a bivariate signal represented as z(t) = $t{x + (t)} + i^t{y + (t)}. In terms of 
the parameters of the modulated elliptical signal of Section |IIlJ the amplitudes and phases of the Cartesian analytic signals are 

a x (t) = K (t)y/l + \\(t)\cos20(t) (70) 
Oy(t) = K{t)y/l - |A(t) | COS 26>(t) (71) 

+3 (In {a(i) cos (9(f) + ir z \b(t)\ sia0(t)}) (72) 
0y(t) = 4>(t)-r z %/2 

(ln{|6(i)|cos0(*) + ir z a(t) sin0(i)}) (73) 

which simplifies the presentation of [12]. Note that the combination 3 (In {z}) is used to implement the four-quadrant inverse 
tangent function, with tan [3 (In {z})] = 9(z)/5ft(z). As \X(t)\ approaches zero, and the signal is nearly circular, both the 
Cartesian amplitudes a x (t) and a y (t) approach K,(t), Meanwhile the phases <p x (t) and [4> y (t) + r z n/2] both approach <fr(t) + 
r z 6(t), which has been previously identified as the phase of whichever rotary component, z + {t) or has the larger 

amplitude. 
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Differentiating the phase expressions (|72| - |73T l one obtains the following, rather complicated, forms for the Cartesian instan- 
taneous frequencies 

,2/ 



u> x (t) 



Uy(t) 



d(j> x {t) _ 
dt 

LUg{t)^l-V{t) 

_ d(j>y{t) 



1 



sin 20(f) d\X(t)\ 



dt 



2^1 

Uj>(t) + r z 



A»(t) 
K 2 (f) 



dt 



(74) 



We (f)^l-A 2 (f) + - 



sin 20(f) d|A(f)| 



(/-/- 



To derive these, we first note that the derivatives of ( |72H73| > become 

d ' 1 



uj x (t) = w$(t) + 3 j— In 
w v (t) = w^f ) + 3 j J- In 



1 



iMtan0(f) 
o(f) 

tan 0(f) 

&(*) 



after pulling out the real parts of the quantities inside the natural logarithms. Then using 

1 d|A(f)| « 



dt 



■In 





dt 




-|A(f)| 








-|A(f)| 


1 



1-A 2 (f) df 



dt 



In 



a(f) 



6(f) 



(75) 

(76) 
(77) 

(78) 



we find ( |74Tl75T l follow in a few straightforward lines of algebra. 

The Cartesian frequencies have three terms: the orbital frequency u>$(t), a weighted version of the precession rate ujg(t), 
and an orientation-dependent term involving the deformation rate. The last term makes an exactly opposite contribution to the 
two power- weighted instantaneous frequencies a 2 (t)uj x (t) and a^{t)oj y (t). Owing to cancelation of these orientation-dependent 



terms, we readily obtain ( 1621 for the bivariate instantaneous frequency using the definition d49b applied to the Cartesian analytic 
signal vector [x+(f) y+{t)] ■ 

Now for the instantaneous bandwidths, differentiating the logarithm of ( |70l - |7T1 > gives 



v x {t) 

Vy(t) 



dln«(t) 1 K 2 (f) d 



dt 
d\n n(t) 
dt 



2a 2 (f) dt 
1 K 2 (f) d 
2a 2 (f) dt 



[|A(f)|cos20(f)] 
[|A(f)|cos20(f)] 



(79) 
(80) 



as expressions for Cartesian bandwidths v x {t) = 4rlna x (t) and v y (t) = 4flna y (t), respectively. In each of these there are 
three terms, due to amplitude modulation, precession, and deformation, respectively; the contributions from each of the latter 
two terms are orientation-dependent. These expressions together with the multivariate instantaneous bandwidth definition (l53l l 
give d63l l for the form of the bivariate bandwidth. 



Appendix B 
Numerical Method 

This appendix describes the means of generating the estimated elliptical signal z(t ) shown in Fig. QJ) and Fig. [4] The reader 
is asked to refer to [12] for details; here we only give a brief overview together with parameter settings. The wavelet transform 
using a generalized Morse wavelet [35], [36] with parameter choices (3 = 3 and 7 = 3 is performed on two real value time 
series x(t) and y(t). This transform is performed at fifty logarithmically spaced frequency bands ranging from one cycle per 
53 days to one cycle per 2.6 days. 

Wavelet ridge analysis consists of locating a maximum of the transform modulus across scale at each time, then connecting 
these points across time into a continuous chain [10], [11], [34]. The wavelet transform along this so-called "ridge curve" 
constitutes an estimate of a (univariate) modulated oscillatory signal. It was shown by [12] that way the wavelet ridge curves of 
the x(t) and y(t) time series, constructed separately, can be combined to form estimates of a modulated elliptical, or bivariate 
oscillatory, signal. Application of this method to the time series shown in Fig. gives the estimated elliptical signal z(t). 
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